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The dynamics of the nuclear-spin quantum computer with large number 
[L = 1000) of qubits is considered using a perturbation approach, based on 
approximate diagonalization of exponentially large sparse matrices. Small 
parameters are introduced and used to compute the error in implementation 
of entanglement between remote qubits, by applying a sequence of resonant 
radio-frequency pulses. The results of the perturbation theory are tested using 
exact numerical solutions for small number of qubits. 

I. INTRODUCTION 

Different solid-state quantum systems are considered now as candidates for quantum bits 
(qubits). They include: nuclear spins 0-gl, electron spins [f|rf|, electron in a quantum dot 
0, and Josephson junctions HH. For the most effective quantum information processing 
the quantum computer should operate with large number of qubits, L. For understanding 
how quantum computer works and for optimizing parameters of quantum computation it 
is necessary to model and simulate quantum logic operations and fragments of quantum 
algorizms. For numerical simulation of quantum computer dynamics one must solve a large 
set of 2 L linear differential equations on a long enough time-interval, or to diagonalize many 
large matrices of the size 2 L x 2 L . Hence, it is important to develop a consistent perturbation 
theory for quantum computation which allows one to predict the probability of correct 
implementation of quantum logic operations involving large number of qubits, in the real 
physical systems. 



In this paper we describe the procedure which allows one to estimate the errors in imple- 
mentation of quantum logic operations in the one- dimensional solid-state system of nuclear 
spins. Our approach does not require exact solution of quantum dynamical equations and 
direct diagonalization of large matrices. We assume that our quantum computer operates 
at zero temperature, and the error is generated only as a result of internal decoherence 
(non- resonant processes). Our approach provides the tool for choosing optimal parameters 
for operation of the scalable quantum computer with large number of qubits. 



II. DYNAMICS OF THE SPIN CHAIN QUANTUM COMPUTER 



Application of Ising spin systems for quantum computations was first suggested in 
Ref. ||10|| . Today, similar systems are used in liquid nuclear magnetic resonance (NMR) 
quantum computation with small number of qubits [II]. The register (a ID chain of N 
identical nuclear spins) is placed in a magnetic field, 
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where < t < t {n+1 \ n = 1,...,M, B z (z) is a slightly non-uniform magnetic field oriented 
in the positive z-direction, b±\ and are, respectively, the amplitude, the frequency, 
and the initial phase of the circular polarized (in the x — y plane) magnetic field. This 
magnetic field has the form of rectangular pulses of the length (time duration) r^ n ' = t^- n+1 ' — 
t^ n \ The total number of pulses which is required to perform a given quantum computation 
(protocol) is M. The chain of spins makes an angle 9, where cos(6 l ) = l/y/3, with the 
direction of the magnetic field B z {z) to suppress the dipole interaction between the spins 
T2\. 

The Hamiltonian of the spin chain in the magnetic field is, 
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where the index k labels the spins in the chain, J is the Ising interaction constant, Q n = 'yb^ 
is the precession (Rabi) frequency, 7 is the gyromagnetic ratio. The function 0W(t) equals 
1 only during the nth pulse and equals to zero otherwise. 

In order to remove the time-dependence from the Hamiltonian (|2|) we present the wave 
function, *&(t), in the time- interval of the nth pulse, in the laboratory system of coordinates 
in the form, 
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where is the wave function in a reference frame rotating with the frequency, v^ n \ 
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the position |1) and of = 1 if the fcth spin is in the position |0), \p) is the eigenfunction of 
the Hamiltonian H . Below we take ip( n > = = for all n. 

The dynamics during the nth pulse is described by the following Schrodinger equation 
for the coefficients A p (t), 

iA p {t) = (E p - Xp n) )A p (t) - § E M*)> ( 4 ) 
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where we put, % = 1, for the Planck constant, and the sum is taken over the states \p') 
connected by a single-spin transition with the state \p), E p is the eigenvalue of the Hamilto- 
nian Hq. In the representation (Q) each spin flip is accompanied by a change of the phase of 
the wave function by the value ±v^ n h, which compensates the time-dependent phase in the 
perturbation, V^(t), in the Hamiltonian (0). In this case, the dynamics of the coefficients, 
A p (t), is governed by the effective time-independent Hamiltonian, H,( n \ and Eq. (|]) can be 
written in the form, iA p (t) = 7i p n p ) A p /(t). 

The dynamics of the coefficients, A p (t), generated by the Hamiltonian, TL^ n \ can be 
computed using the eigenfunctions, A^ n \ and the eigenvalues, e™, of this Hamiltonian by, 

A m (t n ) = E^o(*n-i)E^ ( n) ^ (n) exp(-<r n ). (5) 

m q 

Since the pulses of the protocol are different, we should operate in the laboratory frame and 
make the transformation of the wave function to the rotating frame before each pulse, and 



the transformation to the laboratory frame after each pulse. If we write the wave function 
in the laboratory frame in the form, 

*(t) = E c i»(*)b> ex pH^)' ( 6 ) 
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then the coefficients, C p (t), in the laboratory frame are expressed through the coefficients, 
Ap(t), in the rotating frame as, 

C p (t) = exp(tS p \ n k)A p (t). (7) 

Here 8^ = E p — are the diagonal elements of the Hamiltonian matrix Ttp2, t = t n _i 
before the nth pulse and t = t n — i n _i + r n after the nth pulse. 

III. THE TWO-LEVEL APPROXIMATION 

We explain below how selective (resonant) transitions, which realize a quantum logic 
gate, can be implemented in the system described by the Hamiltonian (Q). For this, we 
consider the structure of the effective time-independent Hamiltonian matrix, Ti, pp i, (here 
and below we omit the upper index, n.) All non-zero non-diagonal matrix elements are the 
same and equal to —0/2. At Q <C 8uj the absolute values of the diagonal elements in general 
case are much larger than the absolute values of the off-diagonal elements, and the resonance 
is coded in the structure of the diagonal elements of the Hamiltonian matrix, 7i pp i. 

Suppose that we want to flip the kth spin in the chain. To do this, we choose the 
frequency of the pulse to be equal to the difference, v = E p — E m , between the energies of 
the states which are related to each other by flip of the kth spin. In this case the energy 
separation between the pth and the mth diagonal elements of the matrix, Ti. pp i, related by 
the flip of the resonant kth spin is much less than the energy separation between the pth 
diagonal elements and diagonal elements related to other states, which differ from the state 
\p) by a flip of a non-resonant k'th (k' ^ k) spin. In this situation, in some approximation, 
one can neglect the interaction of the pth state with all states except for the state, \m), 



4 



and the Hamiltonian matrix, 7i pp >, breaks up into 2 L /2, approximately independent 2x2 
matrices, 
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where £ p = £ m + A pm , |A pm | ~ J or |A ? 
depends on the positions of (k — l)th and (k + 1) spins, V = —fl/2. 

Suppose, for example, that L = 5 and the third spin [k = 3) has resonant (|A pm | = 0) or 
near-resonant (|A pm | ~ J) frequency (we start enumeration from the right spin with k — 0). 
Then, the block 2x2 will be organized, for example, by the following states: 1 01010) and 
|00010); |01111) and |00111); |00001) and |01001), and so on. In order to find the state, 
\m), which forms a 2 x 2 block with a definite state \p), one should flip the resonant spin 
of the state \p). In other words, positions of iV — 1 (non- resonant) spins of these states are 
equivalent, while position of the resonant spin is different. This approximation can be called 
the two-level approximation, since in this case we have 2 L /2 independent two-level systems. 

We now obtain the solution in the two-level approximation. The dynamics is given by 
Eq. (|^). Since we deal only with a single 2x2 block of the matrix Ti. pp ', (but not with the 
whole matrix), the dynamics in this approximation is generated only by the eigenstates of 
one block. 

The eigenvalues ef\ Gq , and the eigenfunctions of the 2x2 matrix (|8D are (we put 
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where A = y/il 2 + A 2 . Suppose that before the pulse the system is in the state |m), i.e. the 
conditions, 

C m (to) = 1, C p (to) = 0, 



are satisfied. After the transformation (|7|) to the rotating frame we obtain, 



A m (t ) = exj>(-iS m t )C m (t ) = exp(-i£ m t ), A p (t ) = 0. 



The dynamics is given by Eq. ([5]), which in our case takes the form: 
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where t = to + r. Applying the back transformation 
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C m (t) — exp(i£ m t)A m (t), 



and taking the real and imaginary parts of the expression in the curl brackets we obtain, 



C m (to + r) = [cos(Ar/2) + z(A/A) sin(Ar/2)] exp(-irA/2). 
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For the amplitude, A p (t), we have, 

A p (t) = A m (t ) [A^A^ exp (-< e f) T ) + Al^A^ • 

& exp {-i [£ m t - (A/2)r]} sin(Ar/2). 
A 

Applying the back transformation, 
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C p (t)=ex V [i(8 m + A)t}A p (t), 



we obtain, 



C p (t + r) = i(n/A) sin(Ar/2) exp(zt A + ir A/2), 



(12) 



When A = (the resonance case) and Ar = 7r (7r-pulse) Eqs. (0) and (12) describe the 
complete transition from the state |m) to the state \p). In the near-resonance case, when 
A / the transition probability is (here we again put the index n indicating the pulse 
number), 



6 



e n = (fi n /A n ) 2 sin 2 (A n r ri ,/2). (13) 

One can suppress the near- resonant transitions and to make the probability, e n , equal to 
zero by choosing the Rabi frequency in the form (see 2irk- method in Ref. ||12||), 



= |A n |/ViF^i. (14) 

The solutions (|ll]) and (|12|) can also be derived without transformation to the rotating 
frame (see Ref. |12[]). However, as will be shown below, our description allows us to introduce 
small parameters and to build a consistent perturbation theory. Using our perturbation ap- 
proach we will compute the dynamics up to different orders in small parameters, and will test 
our approximate results using exact numerical solution for small number of qubits. Below we 
apply the perturbation theory for analysis of the quantum dynamics during implementation 
of a simple quantum logic gate in the spin chain with large number (L = 1000) of qubits. 

IV. PROTOCOL FOR CREATION AN ENTANGLED STATE 
BETWEEN REMOTE QUBITS 

Here we schematically describe the protocol (the sequence of pulses) which allows one to 
create the entangled state for remote qubits in the system described by the Hamiltonian ([|). 
The initial state of the system is the ground state, 1 00 . . . 00). The first pulse in our protocol, 
described by the unitary operator U\, creates the superposition of the states 1 00 . . . 00) and 
1 10 . . . 00) from the ground state, 

£/i|00...00) = - 7 =(|00...00)+i|10...00)). (15) 
v2 

Other pulses create from this state the entangled state for remote qubits. This procedure is 
described by the unitary operator, U', 

C/'4=(|00...00) + i|10...00» = -^(e^ 1 (|00...00) + e^ 2 110 . . . 01)), (16) 
v 2 v 2 



where <pi and if2 are known phases JT2[ 



Now, we describe the procedure for realization of the operator U = U'U\ by applying a 
sequence of resonant pulses. Each pulse is described by the corresponding unitary operator, 
U n , where n — 1, 2, . . . , 2L — 2. (The total number of pulses in our protocol is M = 2L — 2.) 
The unitary operator of the whole protocol is a product of the unitary operators of the 
individual pulses, U = U2L-2U2L-4 ■ ■ ■ U 2 Ui. The first pulse, described by the operator Lq, 
is resonant to the transition between the states |00 . . . 00) and 1 10 . . . 00) . If we choose 
the duration of this pulse as T\ = tt/(20i) (a 7r/2-pulse) then from Eqs. ( |iT1) and (|12"D 
we obtain Eq. (|15|) . In order to obtain the second term in the right-hand side of Eq. 
(|l6l), we choose a sequence of resonant 7i~-pulses which transforms the state |10. . .00) to 
the state |10...01) by the following scheme: (1000...0) |1100...0) 1 1110 ... 0) 
|1010...0) -> |1011...0) -> |1001...0) -►...-> |100...11) -> 1 100... 01). The frequencies 
of pulses which realize this protocol are: = UJ2L-2, ^/( • 3 ' , = ^2L-s, ^ = ^2L-2 — 2J, 
i/( 5 ) = cj 2 L-4, . . . , i/( 2L ~ 3 ) = too — J, i>( 2L ~ 2 } = oj\. If we apply the same protocol to the 
ground state, then with large probability the system will remain in this state because all 
transitions are non-resonant to the ground state. 

Since the values of the detuning are the same for all pulses, A n = A = 2 J (except for the 
fourth pulse, where A 4 = 4J), in our calculations we take the values of Q n to be the same, 
Q n = Q (n ^ 4), and ^4 = 2Q. In this case, the probabilities of excitation of the ground 
state (near-resonant transitions), e n , are independent of n: e n = e, since e n depends only on 
the ratio |A n /f2 n |. (Here e and Q are the numerical parameter used in simulations presented 
below.) One can minimize the probability of the near-resonant transitions choosing e = 0. 



V. ERRORS IN CREATION OF AN ENTANGLED STATE 
FOR REMOTE QUBITS 

Our matrix approach allows us to estimate the error in the logic gate ([16]) caused by flips 
of non-resonant spins (non- resonant transitions). Consider a transition between the states 
|Z) and |/') connected by a flip of the non- resonant k'th spin. The absolute value of the 
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difference between the Zth and I'th diagonal elements of the matrix Tt^S is of the order or 
greater than 5u, because they belong to different 2x2 blocks. Since the absolute values of 
the matrix elements which connect the different blocks are small, |V| <C Sou, one can write, 

*. = *. m +E , - B r a V? ) ' < 17 > 

q i e q e q , 

where prime in the sum means that the term with q' = q is omitted, ip q is the eigenfunction 
of the Hamiltonian TC, the gth eigenstate is related to the Zth diagonal element and the 
g'th eigenstate is related to the Z'th diagonal element, v qq i = 2V(i/j^ > ' ) \I%,\ijj^) is the matrix 
element for the transition between the states ipj® and the sum over q' takes into 

consideration all possible one-spin-flip non-resonant transitions from the state Because 
the matrix Ti is divided into 2 N ~ l relatively independent 2x2 blocks, the energy, (e^? ; ), 
and the wave function, ip^ {ip^), in Eq. ( PD are, respectively, the eigenvalue and the 
eigenfunction with the amplitudes given by Eqs. (|9|) and (|T0| ) of the effective Hamiltonian, 
Ti, in which all elements are equal to zero except the elements related to a single 2x2 block. 

The probability of non-resonant transition from the state \l) to the state connected 
by a flip of the non-resonant k'th spin is, 

Pw = \(l'\%)\ 2 . (18) 



Only one term in the sum in Eq. (|T7| ) contributes to the probability Pw. When the block 
(§) is related to the near-resonant transition (A ~ J), then from Eqs. (H) and ( |10|) the 
eigenfunctions of this block are, 

$f « [1 - (n 2 /32J 2 )]\l>) + (0/47)1*0 « 10. 

w -(0/4 J) |0 + [1 - (0 2 /32J 2 )]|m') w \m'). (19) 

On the other hand, if this block is related to the resonant transition (A = 0), we have 

*f = (1/V2)(|0 + \m% < = (1/V2W) - \m')). (20) 



In both cases v qq i ~ V, so that 



where we put ~ Si, e q °) ~ Sy\ \k — k'\ is the "distance" from the non-resonant fc'th spin 
to the resonant fcth spin. 

The total probability, Hn-i (here the subscript of \i stands for the number of the resonant 
spin), of generation of all unwanted states by the first 7r/2 pulse in the result of non-resonant 
transitions is, 

N ~ 2 1 / tt \ 2 

^ = " £ |iv-i-fe? ' " = l^J ■ (22) 

The probability of error after applying 2A — 2 pulses is, 

2L-2 ^ 2L-2 

P=l ~0 II (!-^)-o II (l-^n-£n), (23) 
Z n=l Z n=l 

where £i = (operation flT5|)) and the last two terms in the right-hand side of Eq. ([23|) are 
related to the last two terms in the right-hand side of Eq. (ITO) . 



VI. IMPROVED PERTURBATION THEORY 



In the consideration presented above, we used the approximate solutions (19) and (20) 



for the wave functions (§) and (10). A more advanced approach, which also does not require 
a diagonalization of large matrices, we use the explicit forms (||) and (|T0D to express the 
wave functions, iffty and ipg) , of the 2x2 blocks in Eq. (|T7D , 

^ = Atr ] w) + 4(y), = aZ^w) + Af^n. (24) 

Then, we put the functions, ip^' and ipg) , into Eq. (|TTD, and obtain the coefficients, A\, for 
the wave function ip q , 

^ q = Y.Al\m), (25) 

m 

where the sum in the right-hand side contains 2L terms. Using the functions, A^, we solve 
the dynamical equations (|5|) with the energies, e q , computed up to the second order of our 
perturbation theory, 
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= 4 0) + E't^. (26) 

q' e Q ~ e q' 

where the prime in the sum means that the term with q = q' is omitted, e^ ' is defined by 
Eqs. © or 

We call the described above approach the "improved" perturbation theory to indicate 
the difference from the approach considered in the previous section. In this approximation 
each eigenf unction, ip q , of the Hamiltonian, TC, is expanded over 2L (see Eq. (^5|)) basis 
functions, \m), with all other coefficients, A q m ,, being equal to zero. Here we use all possible 
transitions between different 2x2 blocks which include the two-spin-flip transitions: a flip 
of the resonant spin and a flip of a non-resonant spin. The number of non-zero coefficients in 
this approximation is 2 L x 2L. It still can be large for large L and can require large computer 
memory for simulation. As will be shown below, under the condition, Q <C J < Sou, this 
approach (the "improved" perturbation theory) gives the results which practically coincide 
with the exact solution. 

VII. NUMERICAL RESULTS 

All frequencies in this section are measured in units of the Ising interaction constant, J. 
Assume that we are able to correct the errors with the probability less than P = 10~ 5 . Our 
perturbation theory allows us to calculate the region of parameters for which the probability 
of error, P, in realization of the logic gate fll6l) is less than P . In Figs. 1 (a) and (b) we plot 
the diagrams obtained by solution of Eq. (|23D and using the improved perturbation theory. 
In the hatched areas the probability of generation of unwanted states is less than P . One 
can see that two approaches yield similar results. They become practically identical at large 
values of the distance between the neighboring qubits, 5oj. 

In almost all quantum algorithms the phase of the wave function is important. We 
numerically compared the phase of the wave function on the boundaries of the hatched 
regions in Figs. l(a,b) with the phase in the centers of these regions, where Q satisfies 
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the 27rfc-method, and the expression for the phase can be obtained analytically (T^|. The 
deviation in phase is only ~ 0.15%. This is much less that the corresponding change in the 
probability, P, of errors (by several orders). 
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FIG. 1. The probability of generation of unwanted states, P, at different values of 5uj and 
$7. In the hatched regions P < Pq (Pq = 10 -5 ). The region delimited by the dashed line is 
obtained using Eq. (|23|). The hatched region delimited by the solid line is obtained using the 



improved perturbation theory described in Section VI. The position of the point A in Q satisfies 
the 2-7rA:-method (14), f2 = £l( k \ where the values of k are indicated in the figures. L = 10. 



We now analyze the probability of errors as a function of 8u. When the value of 5uj is 
large enough, the probability of error (and the widths of the hatched areas in Q) becomes 
practically independent of Sou. This is because at 5u ^> 1 and at e ^> fi the error is provided 
mostly by e, which is independent of 8u. As a consequence, one can, for example, estimate 
the widths of the hatched areas at 5uj ^> 1 taking into account only the near-resonant 
transitions. To do this we put in Eq. (|23|) the value \i n = for all n and obtain 

2L -3 



^ = ^i- 2 n 3 (i-^)) = ^(i-(i-e) 2i - 3 



(27) 



where £„ = £<!, for all n. The positions of the boundaries in Q in Figs. 1 (a,b) can be 
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obtained from the equation Pb = Pq, where Pb is given by Eq. and e is a function of 
tt (see Eq. (^)). 
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FIG. 2. The minimum value of 5uj, 5uj m i n = 5loa (see Figs. 1 (a,b) for L = 10), required to make 
the error in the logic gate ( |i"6| ) below the threshold Pq = 10~ 5 , as a function of number of qubits, 
L. The value of = fiW satisfies the conditions of the 27rfc-method. The values of k are indicated 
in the figure. 

From Figs. 1 (a) and (b) one can see that even when the condition of 27r/c-method is 
satisfied, the error can be large when 5uj (or the gradient of the magnetic field B z (z)) is 
relatively small. Thus, at 5u < 5u>a, where 5uja is the coordinate of the point A in 5u in 
Figs. 1 (a,b), the error is more than Pq. Since the position of the point A in Q satisfies 
the 27rfc-method, this indicates that even in the case when the near-resonant processes are 
suppressed by the 27r/c method (e = 0) the non-resonant transitions can make the error larger 
than the threshold, Pq. We can not suppress entirely the non-resonant transitions, defined 
by the values of \x = (Q/2uj) 2 and L, like we did with the near-resonant transitions. The 
value of Q cannot be decreased considerably because decreasing of Q makes the quantum 
computer very slow, so that the quantum state can be destroyed by decoherence due to 
possible influence of environment. Hence, one can make the value of fi small by increasing 
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8uj. In Fig. 2 we plot the minimum value of 8u = 5uja as a function of the number of qubits, 
L, which was computed using Eq. (f23|). 
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FIG. 3. The exact solution for the probability of generation of unwanted states, P, computed 
using the parameters which correspond to (a) dashed curve AB in Fig. 1 (a) (obtained using 
Eq. (p3|)), (b) lower boundary of the hatched region in Fig. 1 (a) (obtained using the improved 
perturbation theory), (c) dashed curve AB in Fig. 1 (b) (obtained using Eq. (|23|)), (d) lower 
boundary of the hatched region in Fig. 1 (b) (obtained using the improved perturbation theory). 
The dashed lines indicate the solutions obtained using the corresponding perturbation theory for 
the same parameters. 

One can see that Sua becomes large for large L. Thus, for example, for a protons with 
J/ (2tt) ~ 100 Hz (for estimations in this paragraph we use the dimensional units) with the 
distance between the neighboring spins a = 2 nm, the value 5uo / J = 1000 yield the gradient 
of the magnetic field duj/^a cos 9) ~2x 10 6 T/m. From Fig. 2 one can see that this is the 



14 



minimum value of the gradient of the magnetic field for L max « 155 when Q = and 
Lmax ~ 740 when Vl = fi^ 11 ^ required to make the error less than P = 10 5 . At L > L max , at 
a given gradient of the magnetic field, and at O w k — 5 or 11, the error will be always 
larger than Pq. 

In Figs. 3 (a,b) we test our perturbation theory by using the exact numerical solution 
obtained by a diagonalization of 2 L x 2 L matrices and using Eq. (|5[). One can see that there 
is good correspondence with the exact numerical solution for the results obtained using Eq. 
(p3|), and practically exact correspondence for the solution obtained using the improved 
perturbation theory. The similar correspondence can be demonstrated for other parameters 

(5u, n). 

VIII. CONCLUSION 

We developed the perturbation theory which allows us to estimate the errors in imple- 
mentation of the quantum logic gates by the radio-frequency pulses in the solid-state system 
with large number (1000 and more) of qubits. Our perturbation approach correctly describes 
the behavior of the quantum system in the large Hilbert space (the Hilbert space with large 
number of states) and predicts the final quantum state of the system after action of the 
sequence of pulses with different frequencies. This is possible because in the system there 
are small parameters, characterizing the probabilities e, of the near-resonant transitions, 
and probabilities, /i, of the non- resonant transitions, which are small, e < 1 and /i <C 1, 
when the conditions C J C 5w are satisfied. Our approach allows one to control the 
quantum logic operations in the system with large number of qubits and to minimize the 
error caused by the internal decoherence (non-resonant processes). 
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